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ABSTRACT 

Planetary migration is essential to explain the observed mass-period relation for 
exoplanets. Without some stopping mechanism, the tidal, resonant interaction be- 
tween planets and their gaseous disc generally causes the planets to migrate inward so 
efficiently that they plunge into the host star within the gaseous disc lifetime 1-3 
Myrs). We investigate planetary migration by analytically calculating the migration 
rate and time within self-consistently computed, radiatively heated discs around M 
stars in which the effects of dust settling are included. We show that dust settling low- 
ers the disc temperature and raises the gas density in the mid-plane. This inescapable 
evolution of disc structure speeds up type I planetary migration for lower mass bodies 
by up to a factor of about 2. We also examine the effects of dust settling on the gap- 
opening mass and type II migration, and find that the gap-opening mass is reduced 
by a factor of 2 and type II migration becomes slower by a factor of 2. While dust 
settling can somewhat alleviate the problem of planetary migration for more massive 
planets, the more rapid migration of low mass planets and planetary cores requires a 
robust slowing mechanism. 

Key words: accretion, accretion discs - radiative transfer - turbulence - meth- 
ods:analytical - planetary systems:protoplanetary discs 



1 INTRODUCTION 

Planet formation involves the interaction between pro- 
toplanets and gaseous protoplanetary discs. This in- 
teraction drives radial orbital d rifts for protoplanets 
iGoldreich fc Tremaind Il979l . Il980l . hereafter GT80), pos- 
sibly resulting in the diversity of exoplanetary orbital 
properties captured in the observed mass-period relation 
l|Udrv fc Santos! l2007h . In standard disc models, plan- 
ets lose their net angular momentum by resonant in- 
teraction with their natal discs and migr ate inward . 
This is confirmed by analytical calculations (|Wardl 1 19971 : 
iTanaka et alll2002l . hereafte r, TTW02) and numerical simu- 
lations (iNelson et a l. 2000; Kiev et al.ll200ll : lD'Angelo et al.l 
|2002| . l2003l : iBate et al., 20031 in 2D and 3D discs. It is well 
established that the driving force of planetary migration is 
exerted only at Lindblad and corotation resonances (GT80). 
While Lindblad torques excite and propagate density waves 
in discs (Ward 1997i), corotation torques do not. Due to the 
accumulation of angular momentum around the corotation 
region in disc s, without any removal mechan isms such as 
disc viscosity (|Wardlll99ll : iMassetllioOll . \200i ). the corota- 
tion torque saturates (i.e. becomes ineffective). 
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The problem of migration in the standard disc models 
is that, without some sort of slowing mechanism, this tidal 
interaction between a planet and the gaseous disc is so ef- 
ficient that the migration timescale is much shorter than 
the disc lifetime (~ 1 - 3 Myrs). This fast mode of migra- 
tion, known as type I migration, is applicable for low mass 
planets which cannot open up a gap in discs. Based on the 
core accretion scenario, which is the most accepted theory 
for planet formation, cores of gas giants and rocky planets 
are classified as type I migrators. Massive planets such as 
gas giants can open a gap and undergo the so-called type 
II migration, which is abou t two orders of magnitude slower 
than type I migration (e.g. IWardlll997l ). Population synthe- 
sis models confirm that type I migration should be about 
two orders of magnitude slower than expected in order to 
qualitatively reproduce the observed mass-period relation 
(Ida & Lin 2008; Mordasini et al. 2009). One of the greatest 
challenges in planetary formation and migration is to iden- 
tify what physical process(es) succeeds in reproducing the 
diversity of exoplanets and under what conditions planetary 
systems similar to our Solar System are formed. 

This long-standing problem can be resolved by pay- 
ing close attention to the properties of discs. Tidal torque 
strongly depends on the distributions of gas and temper- 
ature in discs and is sensitive to inhomogeneities in these 
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quantities l|Hasegawa fc Pudrit j l2010al . hereafter, HPIO, 
and references herein). Two sources of disc inhomogene- 
ity have been proposed in the Uterature: dead zones and 
ice lines. Dead zones, regions with high density and low 
turbulence (i.e. wh ere the magne t o-rotational instability 
(MR I) is quenched; iGammid Il996l : iMatsumura fc PudritS 
I2OO6I ). drastically affect these distributions and eventu- 
ally provide two of the most robust slowing mechanisms 
(jMatsumura et ahllioOTl . HPIO). Both of them are the con- 
sequence of the difference in turbulence between the active 
and dead zones. One barrier arises by piling up gas at a 
dead zone's outer boundary which is a consequence of the 
time-dependent, viscous evolution of disc s that have strong 
varia t ions in the turbulent intensity (a; IMatsumura et ahl 
I2OO7I . I2OO9I ). The other is the result of a positive tempera- 
ture gradient (HPIO). This radial temperature gradient in 
the dead zone is produced by the back heating from a dusty 
wall which is left in the active region due to the enhanced 
dust settling in the dead zone iHasegawa fc PudritzllioiObl . 
hereafter, Paper I). The variation in turbulence with radial 
direction is the focus in the above treatments. The verti- 
cal variation in turbulence, resulting in layered st ructures: 
the M RI-active surface and the MRI-dead regions (|Gammiel 
1 19961 ) is also important. Combined with the layered struc- 
tures, the ice lin es act as a barrier for type I migration (e.g. 
Ilda fc LinlbOQSl ). This barrier arises because at the ice lines, 
the ice-coated dust density suddenly increases due to the 
low disc temperature. Consequently, the electron number 
density in the surface drops sharply due to the absorption 
by such dust. Since electrons are argents for gas to be MRI- 
active, the gas density of the MRI-active surface region falls 
and gas piles up in the MRI-dead region. As a result, a ra- 
dial positive gradient for the surface de nsity appears, which 
acts as a barrier for type I migra tion (|Masset et al.|[2006l : 
iPaardekooper fc Papaloizoull2009al '). 

In this paper, we step back and investigate the general 
effects of dust settling on planetary migration in homoge- 
neous discs. Although dust settling is observationally con- 
firmed to be ubiquitous in discs around young stars (Paper I, 
references herein), there is no comprehensive study so far on 
the effect that this has on planetary migration. Furthermore, 
dust settling is one of the importa nt processes for planet for- 
mation (e.g. iBirnstiel et al.|[2010l ). In order to proceed, we 
compute detailed, radiatively heated disc models around an 
M star as described in Paper I, which include dust settling 
and the gravitational force of planets embedded in the disc. 
We have two main reasons to focus on M star systems. First, 
low-mass planets such as Super-Earths are the current pre- 
liminary targets for the ongoing and future observational 
missions such as Kepler. M star systems have the highest 
probability of detection using the transit method. In addi- 
tion, some planets around M stars have already been found 
in or near the habitable zone. Second, Monte Carlo calcu- 
lations become computationally expensive for more massive 
discs (such as those around classical T Tauri stars), making 
M star systems an ideal target. 

We generalize our disc models by using minimum mass 
solar nebula (MMSN) models in which the surface density 
behaves as S oc r "'^''^, as well as models with E oc 
(jScholz et al.ll2007l . hereafter S07) used in Paper I and HPIO. 
We then compare the type I migration of various low mass 
planets (2, 5, 10 M^) within our models. Thus, our paper 



provides a detailed exploration of the effect of the thermal 
structure of discs, as a result of dust settling, on planetary 
migration. 

Our plan of this paper is the following. In § [2l we de- 
scribe and compute our disc models and summarize the im- 
portant points for planetary migration for two disc mod- 
els with different surface density profiles. We use these disc 
models as a background disc structure for analytical torque 
calculations. In § [31 we briefiy describe the Lindblad torque 
and discuss the corotation torque. In §[4] and §(5] we present 
our results for the cases of well mixed and dust settling, re- 
spectively. We discuss the effects of individual and combined 
components such as dust settling and the gravitational force 
of planets. The effects of the surface density profiles on the 
migration are also discussed. In § [S] we discuss other heat 
sources and corotation torques, which are both neglected in 
this paper. Finally, we discuss the effects of dust settling on 
type II migration. In § [T] we present our conclusions. 



2 DISC STRUCTURES WITH DUST 
SETTLING 

We adopt the disc models developed in Paper I. In that 
paper, we computed the density and thermal structures of 
radiatively heated discs around an M star, by solving the ra- 
diative transfer equation by mean of a Monte Carlo method 
in 2D axisymmetric {r-z) discs. The main heat source is 
stellar irradiation that is absorbed mainly by dust in discs. 
In order to make our disc models realistic, the gravitational 
force of planets and dust settling are included. We describe 
the treatment of the gravity of a planet in the next subsec- 
tion. In this paper, we adopt two disc models, one of which 
is an MMSN model (E oc r~^/^), and the other of which is 
an S07 model (E oc r^^). These two column density models 
bracket the likely range of models for observed protoplane- 
tary discs. The primary mass is 0.1 Mq and the disc mass is 
4.5 X 1O~^M0, which is increased by a factor of 10 compared 
with the typical disc mass around M st ars in order to allow 
more massive planets to be formed fe.g. lAlibert et al]|2005l ). 
For the other disc parameters, we refer the reader to Paper 
I. 

Note that our r-z disc configurations are very useful for 
calculations of Lindblad torque which is well represented by 
analytical formulae while they are not for corotation which 
needs hydrodynamical simulations in r-cf> discs. However, the 
corotation torque is at least 3 times smaller than the Lind- 
blad torque in our disc models, which is discussed more in 
the next section. Thus, the inability of capturing corotation 
torque in our disc models is unlikely to change our findings 
significantly. 

In this section, we briefly summarise the features of 
our disc models and the resultant temperature and density 
structures of discs. In addition, we discuss the differences 
between the two disc models, and elaborate on the grid sys- 
tems we devised in order to accurately calculate Lindblad 
torque. 

2.1 The gravity of planets 

In this subsection, we briefly summarise the treatment of 
the gravitational force of a planet in our disc models. For 
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complete discussion, we refer the reader to Paper I. We in- 
corporate the gravitational force of a planet into the stan- 
dard vertical hydrostatic equilibrium equation, which may 
be written as 



our grid system by chang ing A_R and n, a nd found AR = 6h 
and n = 6 are sufficient (|Hasegawall2008l 'l. 
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(2.1) 

where ^ = Mp/M,, Mp is the planetary mass, Af, is the 
stellar mass, is the location of the planet from the host 
star, h is the disc scale height, and the normalization con- 
stant Po is chosen so that the densi ty at z = corresponds to 
the un perturbed density (foUowing lJang-Condell fc Sasselovl 
Thus, the dust and gas distributions are perturbed 
by the presence of a planet, and consequently the thermal 
and density structures are affected by the planet. 

We interpolate the density and temperature of discs in 
the region within the Hill rad i us rn ~ rp{Mp/M,Y^^ follow- 
ing [jan^CondeirASasseig^ l|2005h because the hydrostatic 
assumption breaks down within the region. It is well known 
that the horseshoe structure and/or circumplanetary discs 
are located well inside the Hill radius (e.g. iCrida et al...20 09), 
so that we consider this treatment to be very conservative 
since use of the Hill radius overestimates the effect. This 
interpolation is unlikely to be harmful for calculating the 
Lindblad torque because the resonant positions are pushed 
away from their distance from the planets 2/i/3 ~ th- 



2.2 Grid systems 

The Lindblad torque takes its maximum value at the dis- 
tance 2/1/3 from a planet due to the gas pressure effect. 
It th en decreases with increasing distance from the planet 
(Ward ,1997). Thus, it is important to highly resolve the 
thermal structure of the region in the vicinity of planets 
in order to calculate the Lindblad torque accurately. The 
thermal structure of discs, however, is determined globally 
since discs are heated by their central stars. Furthermore, 
the higher the resolution, the longer is the computational 
time. This is because a larger number of photons is required 
in finer grid systems to avoid random noise produced by the 
Monte Carlo methods. Therefore, carefully constructed grid 
systems are required to optimize the computational time. 

We adopt a logarithmic grid system in vertical direc- 
tion. This is sufficient because the inner region in which 
the torques provide the largest contribution is efficiently 
resolved by finer grids. In the radial direction, however, 
more complicated grid systems are required. We combine 
two types of grid systems. Around the planet, we use an 
equally spaced, fine mesh system while for the rest of the 
disc, we use a logarithmically spaced, coarser mesh. The fine 
mesh is inserted into the vicinity of the planets with its size 
Ai? centered at the location of the planet, while the larger 
remaining region is resolved by the coarser mesh. We fix the 
disc size as 50 an and the total number of grids as 180. Also, 
the region Aii is resolved by 65 grids for any case. 

Furthermore, we recalculate the density structure after 
temperature calculations are completed by using n times 
finer systems mentioned above. This treatment enables one 
to reduce the computational time by a factor of 10 or so. We 
performed a convergence study for the Lindblad torque in 



2.3 Well mixed dust 

Dust is important for heating the disc. We first summarise 
the results of our simulations of the thermal and density 
structures of discs with a planet for the case of well mixed 
dust. Fig. [T]shows the dust, gas and total (gas -I- dust) den- 
sity distributions of the MMSN disc models respectively on 
the top to bottom panels. A 10 planet is placed, with- 
out loss of generality, at 6 an in Fig. [T] Also, the zoomed-in 
versions of disc structures around the planet are shown in 
the right column. Note that the disc temperatures shown on 
every panel are the same because they are defined as the 
mass-averaged temperature of dust with various grain sizes 
(Paper I). 

Our numerical simulations confirm that the mid-plane 
region has much lower temperatures than the surface layer 
for both disc models (Fig. [T]). We found, for both disc mod- 
els, that the temperature of the surface is ~ 80 K at 1 an 
while that of the mid-plane is ~ 25 K there. This arises 
because the surface layer is directly heated by the central 
star while the mid-pl ane region is heated only b y the ther- 
mal emission of dust tohian g fc Goldreich|[l997l ). Also, the 
direct heating of discs results in geometrically flared struc- 
tures, which are seen in both the dust and gas density distri- 
butions. Furthermore, the presence of the planet produces 
a low density region above it due to the compression by the 
gravitational force of the planet. Consequently, the region 
above a planet has higher temperatures while the mid-plane 
region, in turn, has lower temperatures. The compression 
results in a small localized peak in a temperature profile in 
the mid- plane (Fig. [2]). 

The difference between the MMSN and S07 disc mod- 
els is that the temperature in the mid-plane region becomes 
lower in the former model, especially at the region around 
the planet and larger disc radii. We found that the temper- 
ature difference between two disc models is about 15 per 
cent in the mid-plane region at 10 an. This is explained by 
the two combined effects arising from a steeper slope pro- 
file. The MMSN disc model has a denser inner region, and is 
therefore more optically thick. In our disc setup, the MMSN 
model is about 4 times denser than the S07 model in the 
mid-plane region at 1 an. As a result, the inner region in the 
MMSN disc model prevents photons from readily penetrat- 
ing the outer region, resulting in lower temperatures there. 
The higher density in the inner region, in turn, results in a 
less dense, outer region even in the mid-plane region, since 
the total disc mass is fixed. The outer region consequently 
becomes optically thin even to the thermal emission of dust, 
resulting in lower temperatures. Since the difference of the 
optical depth is accumulated with increasing disc radius, 
the temperature difference far from the star becomes larger. 
It is interesting that both disc models result in the almost 
identical temperature profiles at the mid-plane, especially at 
smaller disc radii, and are well represented by r~^^^ (Fig- 
ID). 
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2.4 Dust settling 

We showed that dust settling drives the disc structure from 
flared to flatter shapes in Paper I. Fig. [3] shows that the 
same behaviour pertains to the MMSN model as well. Dust 
settling arises from the size distribution of dust grains. In 
general, larger grains cannot be kept aloft even within a re- 
gion with a high level of turbulence (a = 10~^) because 
collisions with gas become inefficient. The increase of larger 
grains in the mid-plane region strongly reduces the probabil- 
ities for smaller grains there to be exposed to and to absorb 
photons emitted by the star. In the surface layer, however, 
the probability for small grains to interact with the photons 
does not change very much since this region is originally 
optically thin. Taking the mass-average of the dust temper- 
atures for all sizes, the surface layer has higher temperatures 
and the mid-plane region has lower temperatures, compared 
with the well mixed case. These lower temperatures in the 
mid-plane, in turn, reduce the scale height of gas, resulting 
in a flatter shape of gas. Table [T] shows the variations of 
the scale height, the temperatures of the surface and mid- 
plane due to dust settling in per cent at 1 au for two disc 
models. Furthermore, geometrically flatter disc shapes are 
favoured because of the temperature of the mid-plane, which 
is well-represented by r~^^'^ (Fig- HI- T his profile is identical 
to an alytical derivations for fiat discs (|Chiang fc GoldreichI 
1 19971 ). 

Dust settling does not change the effects of the gravita- 
tional force of the planet for the two disc models. Thus, the 
gravity of the planet produces higher temperatures above 
it and lower temperatures in the mid-plane region. A sharp 
peak at the co-orbital radii of the planet is enhanced by the 
combined effects of dust settling and the planet for both disc 
models (Fig. |4|. 

The difference arising from the two power law be- 
haviours of the surface density (s = —1 vs -3/2) is qualita- 
tively similar to that found for the well-mixed dust case (see 
table[l]). Since dust settling makes discs fiatter and increases 
the optical depth of the mid-plane region considerably, the 
difference arising from different slopes of the surface density 
is diminished, compared with the well mixed case. We found 
that the temperature difference between two disc models re- 
duces from 15 to 6 per cent in the mid-plane region at 1 au, 
when dust settling is taken into account. Thus, dust settling 
reduces the interaction with photons emitted from the star 
in the mid-plane region for both disc models while it keeps 
the optical depth thick enough for the thermal emission of 
dust even for the MMSN models. 



3 TIDAL TORQUE 

We focus on type 1 migration because it endangers the sur- 
vival of rocky planets and cores of gas giants in the standard 
disc models. We discuss type II migration in § 16.41 We as- 
sume the mass of planets to be small enough to perform type 
1 migration, but to be large enough to disturb the surround- 
ing disc. We neglect any other back-reaction from planets 
onto the disc structure. Thus, we can reasonably estimate 
the migration time analytically for any given disc structure. 
In this paper, we assume planets to be in circular motion 
in order to simply estimate the migration time. In addition. 




Figure 2. Well mixed dust models. The temperature structure 
as a function of disc radius for a MMSN disc model with a 10 
mass planet at 6 au. The solid line denotes the temperature of our 
calculations. The dashed line denotes the best fit to them. A sharp 
peak is produced at the location of the planet by the compression 
arising from the gravitational force of the planet, indicated by 
an arrow. Also, the compression results in lower temperatures 
around it. 




Figure 4. Dust settling models. The temperature structure as 
a function of disc radius for a MMSN disc model with a 10 Mgj 
mass planet at 6 au (as Fig. [2]l. The best fit profile is identical 
to an analytical derivation for fiat discs. A sharp peak at the 
location of the planet is enhanced due to the combined eff'ects of 
dust settling and the presence of the planet. 



we exclude corotation torque because it is safely neglected 
in our disc models as discussed below. 



3.1 Torque at the Lindblad resonances 

Since the basic formula of the Lindblad torque is well dis- 
cussed in the literature, we summarise our assumptions and 
modifications. The complete description is presented in Ap- 
pendix [A] We ad opt the analytical formula of th e Lindblad 
torque derived bv |jang-Condell"fc" Sasselovl l|2005l . hereafter, 
JS05). In the formula, the (total) Lindblad torque r^(r-p) is 



assumed to be calculated by 



r^(rp) 



dz 



dr^(rp,2) 
dz 



(3.1) 



where dV^ {rp, z) /dz is the Lindblad torque exerted by a 
layer at the height from the mid-plane z. Since the mid-plane 
region has the largest gas density, the torque arising there 
provides the largest contribution to the total. Furthermore, 
dust settling transforms the vertical structure of discs from 
fiared to fiatter shapes, and consequently it strengthens the 
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Figure 1. Fully mixed dust models. The density and temperature structures of a MMSN disc model with a 10 mass planet 
at 6 au. The density is denoted by the colors. The scale is shown in the color bar in units of [g cm~'^]. The temperature is denoted 
by the contours in the unit of [K]. Top: dust density. Middle: gas density. Bottom: total density (gas + dust). Blow-ups of the 
region near the planet are shown in the right column. The thick solid lines denote the Hill radius of the planet on every panel. 
Note that the disc temperature shown on every panel is identical to each other because it is calculated as the mass-averaged dust 
temperatures. The surface layer has higher temperatures while the mid-plane region has lower temperatures. The density above 
the planet is reduced by its gravitational force. This results in higher temperatures and low densities above it. 



Table 1. The variation of the disc properties due to the dust settling at 1 au 





MMSN disc models 


S07 disc models 


the scale height h 


10 % reduced 


15 % reduced 


The temperature of the surface region 


28 % increased 


28 % increased 


The temperature of the mid-plane region 


20 % reduced 


28 % reduced 



torque by accumulating the gas density in the mid-plane 
region (see §[5)|. We note that only the gas density is used for 
the torque calculations. One might suppose that the torque 
arising from the dust density is comparable to that from the 
gas for the dust settling case, but this is not the case, because 
the dust density is about 2 orders of magnitude smaller than 
that of gas (see Fig. [3)l . 

The above treatment is very useful for mimicking the 
reduction of the tidal torque in 3D discs, which is found by 
TTW02. In the other words, the integration effectively ac- 
counts for the dilution of the gravitational force of planets 
for the vertical direction without needing to rigorously solve 



the 3D Euler equations (|Wardl 19881 : 1 Artvmowic J 1993h . Sim- 
ilar approaches have been adopted by Menou fc GoodmanI 
(2004, hereafter, MG04) in which the gravitational force of 
planets is weakened by the disc scale height h. Thus, this 
formula includes only 2D modes and relies on the assump- 
tion that the vertical modes arising from a purely 3D effect 
are negligible. Th is assumption is also supported because 
iLubow fc Qgilvid (|i998) who investigated the excitation and 
propagation of density waves at the Lindblad resonances in 
thermally stratified discs, and showed that 2D modes carry 
out more than 95 per cent of angular momentum and verti- 
cal modes are not dominant. 
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Figure 3. Dust settling models. The density and temperature structures of a MMSN disc model with a 10 mass planet at 6 
au (as Fig. [TJ. Top: dust density. Middle: gas density. Bottom: total density (gas + dust). Dust settling makes the surface layer 
hotter and the mid-plane region cooler, resulting in a geometrically flatter disc shape. The presence of the planet provides similar 
effects on the density and temperature structures of the disc. 



For practical purpose, we modified the JS05 torque 
formula. In this modification, the torque density (not the 
torque itself) is calculated, so that one c an avoid calcu- 
lating the resonant positions (|Wardl [l997l ). Consequently, 

dV^ {vp, z) /dz is given as 

where d^V^ (r, z) / dzdr is the torque density at r = r and z = 
z. Under this modification, the Lindblad resonant positions 
r are incorporated into the wave number m (see equation 
HA1|| ). Since the resonant position s or m are controlled by 
the gas pressure (|Artvmowicj|l993h . the lower temperatures 
of the mid-plane region induced by dust settling push the 
resonant positions closer to the planet. Consequently, one 
expects the Lindblad torque becomes stronger, which we 
will confirm (see §[5]). 

3.2 Torque at the corotation resonances 

The torque arising from the corotation resonances is trickier 
to compute because of the possibility of saturation (i.e. it 
becomes zero) as well as having an interaction region that is 



characterized by horseshoe orbits. However, our preliminary 
calculations showed that the unsaturated corotation torque 
in our disc models is at least 3 times s maller than the Lind- 
blad torque l|Paardekooper et al.ll2O10l ). In this calculations, 
we adopted the torque formula consisting of the linear Lind- 
blad torque and both the entropy- and vortensity-related 
horseshoe drags (see their equation (47)). We assumed the 
surface density S oc r"^''^ or E oc r^^ and the disc tem- 
perature T oc r~^/^ with 7 = L4. With the reduction of 
the corotation torque due to (partial) saturation, it is rea- 
sonable to consider that migration is controlled only by the 
Lindblad torque (see Appendix [B] for complete discussion). 
For these reasons, we assume that the corotation torque can 
be safely neglected to zeroth order in our calculations. 

3.3 Migration rate & time 

We calculate the migration rate and time analytically by us- 
ing the tidal torques above, following MG04 and JS05. Thus, 
our calculated migration rate and time are instantaneous 
values rather than time-dependent ones. Since we focus on 
type I migration and assume that the back-reaction from a 
planet on the disc due to the total tidal torque is negligible, 
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(3.3) 



the rate of change of angular momentum of the planets is 

dt 

where J — Mp^GM,rp is the angular momentum of the 
planet. Note that a minus sign appears since the planets 
lose angular momentum by exerting the tidal torque on 
the discs. 

Since the radial drift is due to the migration, we assume 
r-p — rp{t). The migration rate is written as 

(3.4) 

Tp dt J ^ ' 

Since is generally positive, dvp/dt becomes negative, im- 
plying that planets migrate inward. When F^ becomes neg- 
ative, the migration is outward. 

The migration time r,„ig is estimated by 



^l_^y^A4v^ (3.5, 



Again, the migration time becomes positive with F^ posi- 
tive, resulting in inward migration. 



4 RESULTS FOR THE WELL MIXED DUST 
CASE 

We present our results for the well mixed case for our two 
disc models in this section. The density and thermal struc- 
tures of discs discussed in § [2] are used as background disc 
structures for analytical calculations. 

4.1 Torque 

In this subsection, we focus on the behaviour of the torque, 
dr^ {rp, z)/dz, exerted by a layer of gas at height z (see 
equation (|3.2|l 'l. in order to examine the effects of planetary 
gravity on it. Since we found that the basic features of the 
results for the two disc models are similar, we first discuss 
the results for the case of S07 (E oc r~^) in detail and then 
compare the difference between them. 

Fig. [S] shows the torque dT^/dz as a function of z (dis- 
tance from the mid-plane) for the 807 disc models. Except 
for a small dip, the torque dV^/dz of all planets mono- 
tonically increases toward the mid-plane. This trend arises 
partly because the temperature decreases toward the mid- 
plane and partly because the gravitational force of a planet 
is diluted for the vertical direction. The former reduces the 
gas pressure with decreasing z, causing the resonances to 
move closer to the planet. The latter increases the magni- 
tude of the torque with decreasing z. Thus, both physical 
behaviours make the torque of a layer stronger with decreas- 
ing z. The presence of the dip is explained by the distortion 
of disc structure by the presence of the planet, which is dis- 
cussed below. 

As expected, the most massive planet exerts the largest 
tidal torque on the disc (but its mass is still less than the 
threshold mass to open a gap). Consequently, the planet per- 
forms the fastest migration (Fig. [7]). As planetary masses de- 
crease, the tidal torque becomes weaker, resulting in slower 
migration. 

When the MMSN disc model is adopted, the position 



of the dip in the torque moves slightly upward. We found 
that the torque dr^{rp,z)/dz for the MMSN disc model 
is a steeper function of z. In the other words, it becomes 
stronger in the mid-plane region while it becomes weaker in 
the surface layer (also see the left column of Table [2]). This 
occurs because the MMSN model results in lower temper- 
ature structures around the mid-plane region. In addition, 
these lower temperatures cause the gas distribution in the 
disc to be slightly less flared since the gas distribution is de- 
termined by the hydrostatic equilibrium (which is controlled 
by the disc temperature in the mid-plane). 

In order to better understand the presence of the dips, 
we compute the scaled torque exerted by a layer of gas at 
height z, following JS05. JS05 first studied the torque be- 
haviour in radiatively heated discs in which the influence of 
planetary gravity is taken into account. The scaled torque 
is defined as 



dFf 



dz 



1 

1^ 



' dr— {- — 

dz \ p dr 



-JT - (4-1) 



where 



F - i^" 



(4.2) 



and where any quantity with suffix p takes a value at r = rp. 
The scaled torque is designed to isolate the temperature 
dependence of the torque from each layer, since the effects 
of the density distortion on the torque are canceled out by 
the term E/p (see equation (|4.1|l l. 

Fig. [S] shows the scaled torque as a function of z (dis- 
tance from the mid-plane) for the S07 disc models. As with 
the case of the torque dV^ /dz, the scaled torque of all plan- 
ets monotonically increases toward the mid-plane, except for 
the dip. This is consistent with JS05. Due to the Fq in equa- 
tion (|4.1|l , the scaled torques for all planets are normalised to 
the same value at z = 0. It is interesting that the structure of 
the dips is enhanced in the scaled torque. This enhancement 
provides a clear demonstration that dips are located slightly 
below z = th (which is denoted by the vertical lines). This 
is similar to the findings of JS05. We can explain the dip for- 
mation, as follows. The compression of the material above 
the planet due to its gravity, results in higher temperatures 
in the region by reducing the optical depth there. This re- 
sultant temperature structure pushes the resonant positions 
further away from the planet and eventually decreases the 
torque there. Thus, we can conclude that the presence of 
the dip is explained by the distortion of disc structure by 
the presence of the planet. 

Since the temperature distortion which arises from the 
density distortion becomes larger with the mass of the 
planet, more massive planets exert a weaker scaled torque 
on the disc. Thus, the temperature effects arising from the 
distortion diminish the (scaled) torque which is reduced even 
more by massive planets. 

We emphasize that the compression of the material 
above planets, in turn, lowers the temperature of the mid- 
plane region (see §[5|. This makes the torque stronger there. 
Fig. |6] shows that the scaled torques of planets with differ- 
ent masses are almost identical around the mid-plane region. 
This implies that the temperature variation around the mid- 
plane is very insensitive to planetary mass. Since the total 
torque is calculated by summing dV^ /dz from each layer (see 
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Table 2. The difference in the migration time between the S07 
and MMSN disc models at 1 au 





well mixed 


dust settling 


2 Me 


4.0 


1.4 (1.5) 


5 Me 


4.0 


1.4 (1.5) 


lOMe 


3.4 


1.5 (1.8) 



The maximum difference found in discs is given in the brackets if 
the number in the table is not 

equation [3TTJ, the effects of planets on the torque are deter- 
mined by the decrement created by the dips as well as the 
increment around the mid-plane. In our both disc models, 
the increment dominates over the decrement. 

As with the results of the torque dV^ /dz, the main dif- 
ferences between two disc models (MMSN vs S07) are the 
position of the dips and the steepness of the scaled torque 
as a function of z. These differences arises from the change 
in the temperature structure discussed in § 3.2, that is, the 
steeper surface density slope provides lower temperatures 
around the planet and the mid-plane region. Except for 
these, the basic features are the same as above. 

4.2 Migration time &: rate 

In both disc models, massive planets migrate more quickly 
than low mass planets. This trend is clearly shown in Fig. 
[7] which shows the migration time and rate as a function of 
orbital radius for planets with various masses (top and bot- 
tom panels, respectively). The migration time increases for 
planets with large orbital radii, as expected. Also, migration 
in the MMSN disc models is more rapid relative to the S07 
disc models, as discussed above. We found that migration 
in the MMSN disc models is up to 4.0 times faster at 1 au 
than in the S07 disc models (see the left column of Table [2]). 

It is interesting that the migration time for plan- 
ets in discs around M stars is comparable to that of 
classical T Tauri systems. This can be shown by mak- 
ing an order of magnitude est imate. The tidal torque i s 
scaled by Ep{rp/hpf (TTW02: [Paardekooper et al.ll20ld ). 
At 1 au, discs around classical T Tauri systems typi- 
cally have ^„(r„/h„f ~ 10* g cm'^ x (1/0.1)^ = lO'' 
(jChiang GoldreichI [l997l ). For discs around M stars, 
T,p{rp/hpf ~ 10^ g cm-2 x (1/0.01)^ = lO'' (S07). Thus, 
the planetary migration time for discs around both stars is 
similar. Note that, in these above simple calculations, we 
used characteristic values of r-p/hp to estimate an order of 
the timescale. 

4.3 Mass dependence 

In order to further elucidate the effects of the mass of the 
planet, we replot the migration time as a function of plan- 
etary mass at two different orbital radii (Fig. [S]). If these 
effects arising from the gravitational force of a planet were 
neglected, the migration time would scale as an inverse func- 
tion of their masses denoted by the dotted lines (also see 
equations (|3.5|) . (|A4|) and (|A11I) ). Obviously, our calcula- 
tions show some deviation from this function for both disc 
models (although it is not so large). The deviation increases 
with planetary masses, resulting in faster migration. This 



indicates that the increment of the torque around the mid- 
plane is larger than the decrement by the dip formation 
around z = th, and that the difference increases with plan- 
etary mass. 

JS05 found that the deviation increases with planetary 
mass, but it results in slower migration. This indicates that, 
in their disc models, the decrement is larger than the incre- 
ment. It is well known that the dip formation strongly de- 
pends on the grazing angle at whi ch photons emitted from 
the star strike the disc surface (|jang-Condell &: Sasselovl 
l2004h . For the discs around the classical T Tauri stars JS05 
examined, the scale height is large, and consequently the 
grazing angle is also large, compared with the discs around 
M stars we examined. As a result, the decrement by the dip 
formation is enhanced. 

In addition, the deviation becomes smaller with increas- 
ing orbital radius of the planet. This is understood by the 
power-law function of the surface density. Since the disc den- 
sity decreases with the distance from the central star, the 
density and temperature distortions caused by the gravita- 
tional force of the planet are diminished at larger disc radii. 
Consequently, the migration time approaches the simple an- 
alytical solution, which is different from JS05. We suggest 
that the main reason for the difference is that they targeted 
massive discs around classical T Tauri stars including vis- 
cous heating (although our target is low mass discs around 
an M stars without viscous heating). In addition, the orbital 
radii of the planet they considered are 0.5 to 4 au. For such 
a short range, the effects arising from the planet on the disc 
are not diminished. Our results, by contrast, pertain to geo- 
metrically thin, low mass discs and the outer part of massive 
discs for which viscous heating is safely neglected while the 
results of JS05 are for geometrically thick, (the inner part 
of) massive discs in which viscous heating is not negligible. 

It is interesting that the deviation from 1/Mp is almost 
identical for both disc models if the migration times are 
appropriately scaled. Also, the difference of the migration 
times between two disc models becomes small with increas- 
ing the orbital radius of the planets. 



5 RESULTS FOR THE CASE OF DUST 
SETTLING 

In this section, we examine the effects of dust settling in the 
presence of planets with various masses. We also discuss the 
difference between the well mixed and dust settling cases. 

5.1 Torque 

We again focus on the behaviour of the torque, 
dV^ (rp, z) /dz, exerted by a layer of gas at height z (see 
equation H3.2[) ). As is the case of well mixed dust, the re- 
sults of the two disc models are similar, so that we only 
present the results of the the S07 disc model (E cx r~^). We 
discuss the differences between them later. 

Fig. [9] shows the torque dT^{rp, z)/dz as a function of z 
for the planets with various masses. Compared with the well 
mixed case (Fig. [5|, dust settling provides a stronger torque 
around the mid-plane region, as expected (see § This 
arises from lower temperatures in the mid-plane region in- 
duced mainly by dust settling. Since a small change at such 
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z (AU) 

Figure 5. Well mixed dust models. The torque exerted by planets 
on the disc, as a function of the distance from the mid-plane for 
an S07 model (S oc r~^). The solid line denotes the case of 10 
M0, the dotted line is for 5 M0, the dashed line is for 2 M^. 
Every planet is placed at 3 au. Every case shows the dips formed 
by the presence of the planet. The torque is an increasing function 
of planetary mass, as expected. 




z (AU) 

Figure 6. Well mixed dust models. The scaled torque as a 
function of the distance from the mid-plane for an S07 model 
(E oc r~^). The solid line denotes the case of 10 Mgj, the dotted 
line is for 5 Mg, the dashed line is for 2 Mq. Every planet is 
placed at 3 au. The Hill radius for each planet is represented by 
the corresponding vertical line. The scaled torque takes a mini- 
mum value slightly below the Hill radius due to higher tempera- 
tures caused by the gravitational force of the planet. 



lower temperatures strongly affects the resonant positions, 
the distribution of the torque slightly fluctuates, especially 
around the mid-plane region. Although one might worry 
that such fluctuations cause spurious effects, we checked that 
the resultant migration time and rate converge in our grid 
system (Hasegawa 2008). The common features with well 
mixed dust are that the tidal torque becomes stronger with 
increasing planetary mass, and that its distribution shows 
the effects of the planet (although they are small). 

When the MMSN disc models are used, the torque be- 
comes stronger because of lower temperatures around the 
planet, as discussed in the well mixed case (also see the 
right column of Table [21). Thus, dust settling accelerates the 
migration for both disc models by up to a factor of 1.8. (see 
Fig. [11] and Table [3|. 
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Figure 7. Well mixed dust models. The migration time and rate 
as a function of orbital radius of a planet on the top and bottom 
panels, respectively. The solid line denotes the case of 10 Mgj, 
and the dotted line is for 5 Mg), the dashed line is for 2 Mgj. 
The results of the MMSN disc models are denoted by the red 
lines while those of the S07 one are by the black lines. For both 
disc models, the migration time becomes smaller and the rate 
becomes larger as the planetary mass increases and the orbital 
radius decreases. The MMSN disc models produce more rapid 
migration. 




Planetary mass (M^^^tf,) 



Figure 8. Well mixed dust models. The migration time as a 
function of planetary mass. The solid lines with asterisks denote 
the case of rp = 1 au, the dashed lines with open triangles are 
for Tp = 6 au. The results of the MMSN disc models are denoted 
by the red lines while those of the S07 one are by the black lines. 
The two dotted lines denote for a function oc 1/Mp. For both 
disc models, the results deviate from this function due to the 
presence of the planet. The deviation increases with planetary 
masses. Also, our results approach the function as the orbital 
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Table 3. The ratio of the migration time without to with dust 
settling at 1 au 





MMSN disc models 


S07 disc models 


2 Me 


1.3 (1.5) 


1.5 (1.5) 


5 M® 


1.4 (1.5) 


1.4 (1.4) 


lOMe 


1.5 (1.8) 


1.5 (1.7) 



The maximum difference found in discs is given in the brackets 

In order to disentangle the thermal effects on the torque, 
we again consider the scaled torque dV^^^i/dz (see equation 
(14.11) and Fig. I10|) . Compared with the well mixed case, dust 
settling enhances the structure of the dips and moves them 
slightly closer to the mid-plane. This is because dust set- 
tling exaggerates the effects of the planets by transforming 
the disc from flared to flatter structures. Thus, dust set- 
tling diminishes contribution from the scaled torque around 
z = th because of the higher temperatures. At the same 
time, it strengthens the scaled torque around the mid-plane 
region by dust settling as well as planetary gravity. As with 
the well mixed case, the increment of the torque around the 
mid-plane more than compensates the decrement around 
z = vh, and consequently dust settling more accelerates 
planetary migration (Table [3]). 

Compared with the S07 disc model, the MMSN disc 
models tend to move the dips away from the mid-plane 
while the values around the mid-plane region slightly in- 
crease. These trends are the same as the well mixed case 
as already discussed. Compared with the well mixed case, 
however, the scaled torque dVj^^^i/dz for the case of dust 
settling is more sensitive to the change in the background 
disc structure. This occurs because dust settling produces 
cooler and denser disc structures. 

5.2 Migration time & rate 

Fig. [TT] shows the migration time and rate as a function of 
orbital radius of a planet. As discussed above, dust settling 
leads to more rapid migration for both disc models. Table [3] 
summarises the difference of the migration time between the 
fully mixed and dust settling cases for both disc models at 1 
au. Interestingly, the difference between the well mixed and 
dust settling cases is almost identical for both disc models. 

When the results of the two disc models (MMSN vs 
S07) are compared, the difference diminishes (see the right 
column of Table [2]). This arises from that dust settling gives 
similar temperature structures for the two disc models, as 
discussed in § 12.41 




0.0 0.1 0.2 0.3 

z (AU) 



Figure 9. Dust settling models. The torque exerted by a planet 
on the disc, as a function of the distance from the mid-plane for an 
S07 disc model (S oc r~^) (as Fig.fSJl. All three planets are placed 
at 3 au. As is the case of well mixed dust, all three lines show 
the dips formed by the presence of the planet and smoother be- 
haviours. Compared with the case of no dust settling, the torque 
becomes stronger, especially around the mid-plane region and it 
decreases rapidly. 
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Figure 10. Dust settling models. The scaled torque as a function 
of the distance from the mid-plane for an S07 disc model (S <x 
r~^) (as Fig. IS}. All three planets are placed at 3 au. Dust settling 
increases the scaled torque, especially in the mid-plane region 
due to lower temperatures. In addition, it shifts the position of 
a minimum value closer to the mid-plane and enhances the dips 
formed by the presence of the planet, compared with the well 
mixed dust case (see Fig. [gjl. 



5.3 Mass dependence 

In this subsection, we briefly discuss the combined effects 
of planets and dust settling on the torque, because their 
behaviors are the same as those for the well mixed case. We 
confirmed that the deviation from the function (oc 1/AIp) 
becomes larger with increasing planetary mass, and hence 
planetary migration becomes faster. This implies that the 
distortion created by a planet is similar to dust settling even 
though the former is local and the latter is global. Both 
compress the disc density and lower the temperature in the 
mid-plane region. 



Also, we checked that the deviation generally becomes 
smaller with increasing the distance from the central star. 
This arises from the combined effects of dust settling and the 
power-law structure of the surface density. The flatter disc 
shape, (which is a result of dust settling,) reduces the grazing 
angle, resulting in a small decrement by the dip formation. 
At the same time, the power-law structure diminishes the 
effects of planets as discussed at the well mixed case. 
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Figure 11. Dust settling models. The migration time and rate 
as a function of orbital radius of a planet on the top and bot- 
tom panels, respectively (as Fig. [Tjl. For both disc models, the 
migration time becomes smaller and the rate becomes larger due 
to dust settling. As with the case of well mixed dust, the MMSN 
disc models result in more rapid migration. 



6 DISCUSSION 

6.1 Other heat sources 

We have neglected viscous heating of discs in our models. We 
have assumed i nstead that stellar irradiatio n plays the main 
role in heating f Chiang fc Goldreichlll997l l. As discussed in 
Paper I, viscous heating dominates stellar irradiation only 
within 0.1 an for M star systems while the turning point 
is shifted to a few au for cla ssical T Tauri systems (e.g. 
Ijang-Condell fc Sasselov|[200i V Thus, viscous heating can- 
not affect our findings for M star discs beyond 0.1 au. 

In addition, we have neglected the accretion luminos- 
ity of planets. This heat source arises from the accretion 
of surrounding gas, dust, and planetesimals onto the plan- 
ets driving the release of their extra gravitational energy. 
As discussed in Paper I, the resultant temperature profile 
roughly goes to r~^^'^ due to the inverse square law. This 
higher temperature in the vicinity of the planets can slow 
down their migration since the resonant positions are pushed 
away from them due to this higher pressure. In order to ad- 
dress the issue precisely, further detailed numerical studies 
are needed because only a realistic density distribution of 
circumplanetary discs allows one to accurately calculate the 
temperature. Also, such density and temperature structures 



affect corotation torques, which have been neglected in this 
study. 



6.2 Neglect of corotation torque 

Ever since iPaardekooper fc Mellem"al l|2006bl ) numerically 
discovered an entropy-related corotation torque or horse- 
shoe drag, which can be potentially larger than Lindblad 
one, the importance of corotation has been reexamined in 
the literature. However, recent progress has revealed that the 
entropy-related torque is not as important as expected (see 
Appendix |B] for more discussion). One of the reasons is that 
the adiabatic approximation, which establishes the entropy- 
relate d torque, readily breaks do wn in thermal equilibrium 
discs (|Paardekooper et al.l l2010t ) . Consequently, the Lind- 
blad torque likely controls planetary migration over larger 
parameter space than the corotation one. It is one of the 
reasons why corotation torque can be safely neglected in our 
disc models. One caveat is that our disc models cannot cap- 
ture the horseshoe structure which is crucial for the horse- 
shoe drag. Therefore, the proper treatment of this region is 
desired for accurately estimating the corotation torque. In 
addition, our axisymmetric assumption may affect the Lin- 
blad resonances as well particularly through the influence of 
the planet's gravitational field. We do not expect that our 
results are changed significantly. 



6.3 Synthesis of results 

In the above two sections, we examined the effects of plane- 
tary gravity and dust settling on the tidal torque in detail. 
In this subsection, we integrate our results and discuss how 
our results map to the literature. 

We pro ceed by comparing our results with those of 
IWardI |l993) and TTW02 (Fig. [H]). We adopt the results 
for the S07 disc model with a 5 planet. Fig. 1121 shows 
an interesting transition behaviour in the tidal torque for 
radiatively heated discs. For planets with small orbital radii 
(r < 3 au), our results for the well mixed and dust settling 
cases are close to the results derived from 2D discs, while, for 
planets with large orbital radii (r > 3 au), ours are close to 
the ones derived from 3D discs. Obviously, this arises from 
the complex combination of the effects of planets on the 
thermal and density structures of discs, as discussed above. 
Thus, in order to accurately estimate the migration time, it 
is essential to consider the effects of stellar irradiation and 
planetary gravity simultaneously. 

For both cases, the migration time is shorter than the 
disc lifetime, which is a order of 10^ years. Indeed, dust 
settling aggravates the problem of type I migration. There- 
fore, it is clear that robust slowin g mechanisms are needed 
to resolve the problem. HP 10 and iMatsumura et al.l (|2009l ) 
demonstra t ed th at dead zones act as barriers for it while 
llda fc LinI (|2008l ) considered that ice lines in layered struc- 
tures play such a role. 



6.4 Gap-opening mass & type II migration 

We have focused so far on type I migration. For complete- 
ness, we now examine the effects of dust settling on gap- 
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Figure 12. The migration time as a function of the distance from 
the star. We adopt the S07 disc models with a 5 M^pi^^ planet. 
Our results are denoted by the black and red solid lines for the 
well mixed and dust settling cases, respectively. For comparison 
purpose, the results of TTW02 for 3D discs are denoted by the 
black dashed line, those of TTW 2 for 2D discs are by the black 
dotted line, and those of IWardI lll997l ) are by the black dash- 
dotted line. Our results show an transition in the migration time 
from the 2D to 3D calculations. 



opening mass planets which undergo type II migration in 
this section. 

Gap-opening mass is well d iscussed in the literature 
fe.g. iMatsumura fc PudrTt3l2006l . references herein), so that 
we briefly present criteria above which planets can open up 
a gap . For inviscid discs, a gap-opening mass is 



> 3C 



hr. 



(6.1) 



We introduce a coefHcient C into the original condition, 
in order to take into account (some) possible reduction of 
the tidal torque. In many previous studies, the gap-opening 
criterion is derived under the 2D approximation for discs. 
TTW02, however, showed that the tidal torque in 3D is a 
factor of 2-3 smaller than that in 2D. Since the gap-opening 
mass involves the tidal torque, the mass for 3D discs can 
increase by at least a factor of 2 or 3. We therefore take 
C = 3 here. We leave a more detail study on it to a future 
publication. 

For viscous discs, a gap-opening mass is 



(6.2) 



where a is the strength of turbulence. The presence of vis- 
cous torque which plays a counteractive role in gap-opening, 
that is, it closes a gap. 

Fig. [13] (Top) shows the gap-opening mass as a function 
of the distance from the star. In general, the mass derived 
from viscous discs is larger than that from inviscid discs. 
Both criteria for the dust settling case are smaller than those 
of the well mixed case, as expected. This is a consequence 
of having flatter discs which arise from dust settling. We 
note that, even for the dust settling case, all planets we 
have considered (2, 5, 10 A%) satisfy the criterion of type I 
migration (i.e. not open-up a gap) for r = 1 au to r = 10 an. 
We found that dust settling reduces the gap-opening mass 



at 10 au, respectively, for inviscid discs by about a factor of 
2.5, for viscous discs by about a factor of 2. 

In the above analys i s, ther e is no time-dependency. 
iPaardekooper fc Melleinal (|2006al ') simulated the interaction 
of gas and dust in discs in the presence of massive plan- 
ets, using a two-fluid approach. Although they cannot in- 
clude dust settling, since they considered 2D discs, they 
demonstrated that dust-gaps are more readily formed than 
gas-gaps. Furthermore, they found that dust removed from 
the location of planets is accumulated at mean-motion 
resonances. Combined with such simulations, gap opening 
masses for dust settling are likely to be much lower. 

Once a planet o pens up a g ap in the disc, it undergoes 
Type II migration (|Ward|[l997l ). This type of migration is 
controlled by the viscous evolution of the gas in discs. The 
timescale is generally written as 



TII = 



(6.3) 



where ftp is the angular frequency at r = rp (|Teraueml2003l l. 
Thus, type II migration is also affected by dust settling via 
the reduced disc scaleheight (the bottom panel of Fig. I13|) . 
We found that dust settling increases the timescale by about 
a factor of 2 at 10 au. 

As a result, dust settling somewhat alleviates the prob- 
lem of planetary migration in two respects: the reduction of 
the gap-opening mass and the slowing of type II migration. 
In discs with dust settling, low mass planets undergo faster 
type I migration. At the same time, such planets more read- 
ily open-up a gap, and switch to type II migration which is 
generally slower than type I. 



7 CONCLUSIONS 

We have carefully examined Lindblad torques in radiatively 
heated discs. In order to make our disc models realistic, we 
have included dust settling and the gravitational force of 
planets. We have demonstrated that dust settling accelerates 
type I migration by up to a factor of about 2. This arises 
from the lower temperatures of the mid-plane region and the 
resultant flatter shape of discs. Both are the consequence of 
dust settling. 

We have also investigated the mass dependence on the 
torque, by considering planets with various masses (2, 5, 10 
M^). We have found that the gravitational force of planets 
has two effects on the torque. One is the reduction of the 
torque contribution around z = th- The other is to lower 
the disc temperature in the mid-plane which increases the 
torque. In our disc models, the increment is dominant over 
the reduction, so that massive planets undergo even more 
rapid migration. 

We have also considered disc models with two differ- 
ent surface density profiles: E cx (S07) and E oc r~^^^ 
(MMSN). We found that, for the well mixed case, planets for 
the MMSN models migrate faster than those for the S07, by 
up to a factor of 4. For the dust settling case, the difference 
becomes small (about a factor of 2). 

Finally, we have studied the effects of dust settling on 
the gap-opening mass and type II migration. The reduc- 
tion of the scale height by dust settling decreases the gap- 
opening mass and increases the timescale of type II migra- 
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Figure 13. The gap-opening mass and tlie timescale of type II 
migration as a function of the distance (on top and bottom panels, 
respectively). We adopt the S07 disc model without planets. The 
gap-opening masses for viscous discs are showrn with solid lines 
while those for inviscid discs with dashed lines. For both panels, 
the results of the well mixed case are denoted by black lines while 
those of the dust settling case are by red lines. The gap-opening 
masses for the viscous discs are larger than those for the inviscid 
discs, as expected. Also, dust settling reduces the gap-opening 
mass, compared with the well mixed case. For the timescale of 
type II migration, dust settling results in slower migration than 
the well mixed case. 



tion. Thus, dust settling may slightly reduce the problem 
of planetary migration for sufficiently massive planets. Our 
results show that dust settling has a clear and significant 
effect on planetary migration of both low and high mass 
planets. It is clear that an even more effective slowing mech- 
anism is required in order to slow type I migration. Cur- 
rently, tw o types of inhomog e neity in discs could do the job. 
HP 10 and Matsu mura et al.! (2009) demonstrated that dead 
zo nes provide two r obust barriers while ice lines considered 
bv lldafcLini d 20081 ) provide another robust barrier. In our 
next paper, we will incorporate these barriers in one disc 
model, and discuss their roles in producing the mass-period 
relation. 
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APPENDIX A: THE LINDBLAD TORQUE 
FORMULA 

Here, we describe the formulation of the Lindblad torque. 
Since we calculate the torque density, the Lindblad reso- 
nances r are incorporated into the wavenumber m; 



i(r,z) = 



(Al) 



is exerted by the planets on discs. Thus, the positive sign 
of e in equation (|A5|) means that the planets exert the tidal 
torques on discs, resulting in the loss of angular momentum 
from the planets and vice versa. 

The tidal torque reaches its maximum value at the reso- 
nant position which is located at about 2/i/3 from the plan- 
ets. It decreases far from this point because the di stance 
between the planets and the resonances is increased (|Wardl 
Il997h . Furthermore, it is negligible within the region be- 
tween ± 2/1/3 due to the gas pressure. In other words, the 
resonant positions are pushed away from the planets due to 
the effe ct of the azimutha l pressure, - the so-called pressure 
buffer (| Artvmowic j [199^ ) . This relation is clearly seen in 
equation I|A1[) . When ~ Hp, m becomes imaginary, mean- 
ing that there is no m corresponding to such resonant posi- 
tion r. Thus, the dominant contribution to the total torque 
arises from the region with r ~ ± 2h/3. 

This allows one to approximat ely replace c os in equa- 
tion dMl) with 1 - 9^/2 (GT80; IWardI Il997l . JS05). The 
Laplace coefficient 67/2 then becomes 



where f2p = ^/GM^/r^ is the angular velocity of a planet 
at the distance Vp from the host star with its mass M, , 
Cs (oc VT) is the sound speed, T is the disc temperature, the 
angular velocity Q, is 



Q.^{r,z) = Q.Kep ( 1 + ^ 



-3/2 



+ 



_l_9p 

rp dr ' 



(A2) 



f^ifep ~ \/ GM, I r'^ is the Keplerian frequency, p is the gas 
density, p = c^p is the pressure, and the epicyclic frequency 
K is 

1 d , 4„2X 



K (r,z) 



^3 



(A3) 



Thus, m becomes a function of r that represents the reso- 
nant positions. Note that the vertical thickness of discs are 
included in equation (|A2[) . 

For a planet with mass Mp, the Lindblad torque density 
d r jdzdr of a layer at the height from the mid-plane z may 
be written as 



d /dF^ 



dz V dr 



{r,z) 



2p pr n 4 2 
V(l + 4a.^ "^^ ' 



where 



e = 



for the outer resonances, 
for the inner resonances. 



(A4) 



(A5) 



p = Mp/M*, 5 = mcs/rhi, the forcing function is 

TV Qr rf6l/2(ar,C) 



m 



dOr 



+ 2yrTe&r/2K,c) 



and the Laplace coefficient 6" 



/2 



COS mOdO 



1" J a y/l - 2ar cos6I + q2 +(^2 



(A6) 
(A7) 



where = r/rp is the normalised Lindblad resonant po- 
sition and C, = z/r diminishes the 6^" 2 due to the vertical 
thickness of discs. Note that the integration of equation (|A4|) 
in terms of z with C = leads to the famous analytic al for- 
mula e of the Lindblad torque density in 2D discs fsee lWardl 
I1997I . MG04). Also, note that the tidal torque shown above 



?'l/2(ar,C) 



^o(A), 



(A8) 



where Ki is the modified Bessel function of the second kind 
of order i and the argument A is 



A(a,.,0 = m 



(A9) 



Subs tituting in this approx i matio n and using Kq = —Ki 
(e.g. lAbramowitz fc StegunI 1 19721 ) , the forcing function 
becomes 



/Qr 

(AlO) 

Finally, the total torque exerted by the planets on discs 
is calculated by 

F^^r^^f d.^(^(.,.)). (All) 



APPENDIX B: THE PROBLEMS WITH 
COROTATION TORQUE 

Here, we discuss corotation torque and summarise its prob- 
lems. The corotation torque was originally derived by GT80. 
The corotation resonances arise when the orbital frequency 
of the gas is identical to the orbital frequency of the planet 
Qp. The analytical formula is 



where 



dr I dr \B 



B=^^(r^n) 



(Bl) 



2r dr 

is one of the Oort constants or the vorticity of gas, and 



(B2) 



GA4 



bT/2{ar) 



(B3) 



is the mth order Fourier component of the potential of the 
planets. Thus, the sign of the corotation torque is deter- 
mined by the radial gradient of vorticity per unit surface 
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density, called the vortensity. For the Keplerian discs, B oc 
i^Kep- Consequently, FjJ^ vanishes with E cx r~^^^. Unfortu- 
nately, this first derived analytical formula had a discrep- 
ancy with numerical calculations l|Korvcanskv fc PoUackl 
1 19931 ) ■ Furthermore, the mechanism of the exchange of an- 
gular momentum between the planets and the gaseous discs 
and its relation with the the gradient of the vortensity is 
physically not clear. 

The discrepancy between the analytical and numeri- 
cal torques was resolved by TTW02 by replacing (fif^-^ with 
1 0m + Vm\'^, where rjm is the mth Fourier component of the 
enthalpy perturbation. This implies that the gas pressure or 
temperature around the planet is crucial for a proper pre- 
scription of the corotation as well as the Lindblad torques 
since the enthalpy is defined by the pressure divided by den- 
sity. The revised full analytical f ormula of the corot ation 
torque has recently been derived bv lZhang &i Lail |2006l ). We 
call the torque discussed above the linear corotation torque 
sinc e it is derived from the linearized equations. 

IWardI l|l99lf) derived the corotation torque based on a 
fundamentally different approach. This corotation torque is 
often called a horseshoe drag to distinguish with the linear 
one. The horseshoe drag arises from the interaction between 
a planet and the fluid element moving in the vicinity of it. 
When the fluid element encounters the planet, the orbital 
radius of the element is shifted from the inner to the outer 
one relative to the planetary orbit. During this non-linear 
process, the fluid element gains angular momentum from 
the planet and vice versa since angular momentum of the 
disc increases with distance from the central star. When the 
orbital radius of the fluid element is shifted from the outer 
to the inner orbits, the sit uation is the o pposite. As seen in 
numerical simulations (e.g. lMassetll200ll ). these shifts result 
in a horseshoe orbit. That is why it is called the horseshoe 
drag. Furthermore, adiabatic invariance forces the fluid to 
continue its circular motion l|Wardlll99j ), which prohibits 
the excitation of density waves and possibly also produces 
saturation (that is, the horseshoe drag becomes negligible). 

ft is very interesting that the horseshoe drag also de- 
pends on the radial gradient of the vortensity which is ex- 
plained by mapping process of the area occupied by the fluid 
element. Note that there is no reason to identify the linear 
corotation torque with the horseshoe drag since they are de- 
rived from different approaches. The formula of the horse- 
shoe drag is 



(B4) 



where Xs is the half width of the horseshoe region and 

is one of the Oort constants (|Wardlll99ll ). Thus, the horse- 
shoe drag strongly depends on Xs althou gh Xa is treated as a 
parameter in any numerical simulat ions (|Masset et al.|[200^ : 

IPaardekooper fc Papaloizoull2009bl'l . 

Recently, Paardekooper fc Papaloizoul l|2009al 'l have re- 
examined the linear corotation torque and horseshoe drag 
and found that the linear one is only valid for the early 
stage of numerical simulations or for discs with very high 
viscosity (a ~ 0.1). Except for these, the total torque is 
well-represented by the sum of the Lindblad torque and the 
horseshoe drag. With the proper choice of Xs, the horseshoe 



drag is about 2-3 times larger (in magnitude) than the linear 
corotation torque, but it is still smaller than the Lindblad 
one for discs with s < 1, where E cx r^. A more complete 
treatment for the horseshoe drag including the effects of the 
pressure is done bv lCasoli fc Massetl (|2009l ') and confirms the 
same results. 

In the above studies, the temperature which may be 
a function of radius from the star is assumed to be fixed 
during mass transfer. In other words, discs are assumed 
to be barotropic. In this case, the magnitude of the lin- 
ear corotation torque or horseshoe dr ag is always smaller 
than the Lindblad one (|Korvcanskv fc P ollack 1993). There- 
fore, one can safely assume both torques to be negli- 
gible even if it is not saturated. The situation, how- 
ever, can change by relaxing th e isothermal assumption 
d Paardekooper fc Mellema l2006bl ). For discs in which the 
energy conservation is fully considered, the linear coro- 
tation torque and horseshoe drag are determined not 
only by the radial gradient of t he vortensity but a l so by 
the radial gradient of entropy ( Baruteau fc Massetl l2008l : 



IPaardekooper fc Papaloizoull200g ). It is interesting that the 
entropy-related corotation torque, especially horseshoe drag, 
can dominate over the Lindblad one and control migration 
if discs have a large entropy gradient and are considered to 
be adiabatic. In the region with a large (in magnitude), neg- 
ative entropy gradient, planets migrate outwards due to the 
large, positive horseshoe drag. The complete analytical for - 
mulae in 2D discs are derived bv lPaardekooper et al.l (|201(J ). 
ass uming saturation d oes no t occur yet. 

iMasset fc Casolil (|2009i ) reexamined the corotation 
torque in adiabatic discs, taking into account the effects 
of pressure. They showed that the entropy-related torque 
does not arise from the overdense and underdense regions, 
as had been previously thought. Instead, they showed that 
entropy gradients make the horseshoe region considerably 
asymmetric, and this asymmetry exerts the entropy-related 
torque by the excitation of evanescent waves at the horse- 
shoe sparatrices. Since these evanescent waves are excited 
by pressure disturbance that are a result of horseshoe orbits, 
the magnitude of the entropy-related torque depends on the 
p erturbed pressure. This giv es an explanation of the results 
of lBaruteau fc Massetl (|2008l ) that the corotation torque in a 
very cold, adiabatic disc is identical to that in an isothermal 
disc. 

Thus, the corotation torque or horseshoe drag can be 
important only for discs assumed to be adiabatic with high 
temperatures, that is, the radiative cooling timescale in discs 
is longer than the planet's orbital period. Otherwise, the di- 
rection of migration is determined by the Lindblad torque. 
The horseshoe drag, however, has some problems. One is 
the width of the horseshoe region, and the other is sat- 
uration. The former is characterised by the mass of the 
planets, but it is parameterized in numerical simulations 

J iartly because a planet is treated as a point m ass object 
D'Angelo et"alll2003l : IPaardekooper et al.llioTol ') and partly 
because an y simulation cannot infinitely resolve the horse- 
shoe region ("M asset et al.|[200^ : |Pa ardekoop er fc Papaloizoul 
[2009b). For the latter problem, the horseshoe drag can be 
saturated (i.e. negligible) due to adiabatic invariance. 

We emphasise that the problem of saturation is essen- 
tially how the disc viscosity connects the horseshoe region to 
the rest of disc by transferring angular momentum between 
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them (e.g. iMasset fc Casolil [ioiol ). This is crucial for un- 
derstanding hors eshoe drags. As a example , we summarise 
the work done bv lPaardekooper et al.l ()201CI ) who considered 
the saturation of the vortensity-related and entropy-related 
horseshoe drags. For the vortensity-related one, the satura- 
tion i s controlled only b y viscous diffusion (also see iMasseO 
I2OOII . I2OO2I : IWardll2007l ). For the discs with a ~ 10-^ the 
horseshoe drag becomes almost zero while it is unsaturated 
for the disc with a ~ 0.01 (see their fig. 2). For the entropy- 
related one, which is more important for migration, the sat- 
uration is_controUed_bj;_viscou^^ the thermal dif- 
fusi on (Paardekooper fc Papaloizoul 12008?). It is important 
that lPaardekooper et al.l l|2010l ) found the thermal diffusion 
is always much larger than the viscous diffusion in thermal 
equilibrium discs in which viscous heating is balanced with 
radiative losses. Since thermal diffusion defines the validity 
of the adiabatic approximation of discs, this finding strongly 
restricts the applicability of the entropy-related horseshoe 
drag. Indeed, they have confirmed that the entropy-related 
horseshoe drag is reduced by about a factor of 2 due to ther- 
mal diffusion (see their fig. 7). Although their findings are 
potentially changed by the mass of a planet which controls 
saturation (IPaardekooper et al.ll2O10l ). we can conclude that 
the entropy-related horseshoe drag is not so important for 
migration as expected before. This is because high viscosity 
(a ~ 0.01) is required for horseshoe drag to be unsaturated 
while it also boosts thermal diffusion, resulting in the discs 
considered to be (locally) isothermal. Consequently, only the 
vortensity-related horseshoe drag matters. For low viscosity 
(a ~ 10~^), any torque arising at the corotation resonances 
readily suffers from saturation. Therefore, it is reasonable 
that corotation torque is assumed to be negligible. 

The bottom line appears to be that the corotation 
torques or horseshoe drags are unlikely to play a dominant 
role in planetary migration. 

This paper has been typeset from a TJrjX/ I^Tp^X file prepared 
by the author. 



